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ABSTRACT 

In the primordial universe, low mass structures with virial temperatures less than 10 4 K were unable 
to cool by atomic line transitions, leading to a strong suppression of star formation. On the other 
hand, these "minihalos" were highly prone to triggered star formation by interactions from nearby 
galaxy outflows. In Gray & Scannapieco (2010), we explored the impact of nonequilibrium chemistry 
on these interactions. Here we turn our attention to the role of metals, carrying out a series of high- 
resolution three-dimensional adaptive mesh refinement simulations that include both metal cooling 
and a subgrid turbulent mixing model. Despite the presence of an additional coolant, we again we 
find that outflow- minihalo interactions produce a distribution of dense, massive stellar clusters. We 
also find that these clusters are evenly enriched with metals to a final abundance of Z « 10 -2 Z . As 
in our previous simulations, all of these properties suggest that these interactions may have given rise 
to present-day halo globular clusters. 

Subject headings: galaxies: formation - galaxies: high-redshift - star clusters: general - globular clusters: 
general - shock waves - galaxies: abundances 



1. INTRODUCTION 

Turbulent processes are instrumental in understanding 
a wide range of astrophysical observations, including the 
elemental homogeneity in field stars (Reddy et al 2003), 
open clusters (e.g., Friel & Boesgaard 1992; Twarog et 
al 1997; Carraro et al 1998), the Magellanic clouds (Ol- 
szewski et al 1991), dwarf- irregular galaxies (Thuan et al 
1995), and galactic and disk H II regions (e.g., Deharveng 
et al 2000, Henry & Worthey 1999). Turbulent mixing 
is also important for the enrichment of primordial gas 
(Pan et al 2007) and the transition from Population III 
to Population II stars (Scannapieco et al 2003). 

Similarly, turbulence affects the distribution and for- 
mation of molecular species. Because most chemical re- 
actions are strongly temperature-dependent, by simply 
moving material to a region with different physical prop- 
erties such as density, temperature, and UV flux, or by 
creating a local heating event through turbulent dissi- 
pation, many reactions can be greatly enhanced, which 
alters the final abundance of each species (for a review 
see Scalo & Elmegreen 2004). This can be particularly 
important in primordial gas at high-redshift, whose cool- 
ing properties are highly dependent on the mass fraction 
of molecular hydrogen and hydrogen deuteride. Further- 
more, turbulent transport is equally important after re- 
actions occur, as it will alter the physical distribution of 
the newly formed chemical species (e.g, Xie et al 1995). 

In the first paper of this series we explored the effect 
of nonequilibrium chemistry and associated cooling in 
the interaction of a high-redshift galactic outflow and a 
low-mass primoridal cloud (Gray & Scannapieco 2010, 
hereafter Paper I), in the absence of turbulence or met- 
als. A large population of high-redshift gravitationally 
bound clouds with virial temperature below 10 4 K is a 
generic prediction of the cold dark matter (CDM) model 
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(Bromm & Larson 2004 and references within). These 
clouds, which are often called "minihalos," are too small 
to excite the necessary radiative transitions to cool via 
atomic hydrogen and helium lines. Thus further mini- 
halo cooling requires either molecular-line transitions or 
metal/dust cooling. Although some small abundance 
of H2 is expected after recombination and important in 
cooling of larger structures (e.g., Abel et al 2002; Bromm 
et al 2002) the resulting Lyman- Werner photons (11.2- 
13.6 eV) from the stars in these objects (e.g., Haiman et 
al 1997; Ciardi et al 2000; Sokasian et al 2004; O'Shea & 
Norman 2007) would likely dissociate this trace amount. 
Furthermore, even if some trace amount of H2 exists it 
is unlikely that it would impact the structure of these 
objects (Whalen et al 2008a; Ahn et al 2009). 

In Paper I we used high-resolution three-dimensional 
adaptive mesh refinement (AMR) simulations to model 
the interactions between minihalos and high-redshift 
galactic outflows, including a 14 species model of pri- 
mordial chemistry. We found that several dense clusters 
with masses above 10 4 M were formed, which may re- 
semble present day halo globular clusters. 

Here we consider the role of turbulent mixing on these 
interactions, returning to this model and including the 
effects of unresolved turbulence and metal-line cooling. 
In particular, we are interested in three primary ques- 
tions: First, how does the inclusion of turbulence affect 
the mixing and final abundance of metals coming from 
the galactic outflow? Second, does metal-line cooling sig- 
nificantly alter the final state of the cloud? And finally, 
do we obtain a similar distribution of clusters as our pre- 
vious simulations without subgrid turbulence or metals? 

The structure of this paper is as follows. In Section 
2 we describe our approach to the subgrid modeling of 
turbulence, compare the model to analytic test problems 
and laboratory experiments, and describe our method 
for metal line cooling. In Section 3 we outline our model 
setup and initial conditions, and in Section 4 we present 
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our results and their implications to local observations. 
Conclusions are given in Section 5. 

2. NUMERICAL METHOD 

All simulations were performed using FLASH ver- 
sion 3.1.1, a multidimensional adaptive mesh refine- 
ment hydrodynamic code (Fryxell et at 2000) which 
solves the Riemann problem on a Cartesian grid using 
a directionally-split piecewise parabolic method (PPM) 
solver (Colella & Woodward 1984; Colella & Glaz 1985; 
Fryxell et al 1990). 

In Paper I, we developed and verified a 14 species 
chemical network that traced the evolution of both 
atomic (H and He) and molecular (H 2 and HD) species, 
which included all the pertinent cooling rates. In the 
current work, in order to study the evolution and im- 
pact of metals in outflow-minihalo interactions we add 
two further packages to our simulations: a turbulence 
model that tracks the subgrid mixing of enriched and 
primordial gas and a cooling module that accounts for 
additional cooling in the presence of metals. Here we 
describe each of these in turn. 



2.1. K-L Turbulence Model 

Within FLASH we have implemented a buoyancy and 
shear driven model of turbulence using a two equation K- 
L model, where K represents the turbulent kinetic energy 
and L represents the eddy length scale. The model was 
originally developed and used with great success to de- 
scribe turbulent fluid flow in inertial confinement fusion 
(ICF) experiments (Dimonte & Tipton 2006, hereafter 
DT06; Chiravale 2006), and it reproduces three primary 
fluid instabilities: the Ray leigh- Taylor (RT) instability, 
which arises when a low density fluid supports a high 
density fluid under the influence of gravity or acceler- 
ation, the Richtmyer-Meshkov (RM) instability, which 
occurs when a shock interacts with a fluid of different 
acoustic impedance such as a density gradient, and the 
Kelvin-Helmholtz instability, which arises from velocity 
shear between two fluids, even when they have otherwise 
identical properties. 

In Scannapieco & Briiggen (2008; hereafter SB08), the 
authors used the original DT06 model to study AGN- 
driven turbulence in galaxy clusters focusing purely on 
RT and RM driven turbulence. For our expanded model, 
we have added the additional contribution from the KH 
instability, which is crucial to the problem we are study- 
ing. The DT06 model is based on the Navier-Stokes 
equations expanded to include a turbulent viscosity \±t 
and pressure Pt which are dependent on the eddy size L 
and the turbulent kinetic energy K. To compute these 
properties we divide the flow into two components; writ- 
ing velocity, for example, as the sum of mean u and fluc- 
tuating u" components: 



3D, 



u + u" 



(1) 



where u is the mass averaged variable u = ~pu/p, pu" = 
0, p is the mass density, and the overbar represents an 
ensemble average over many realizations of the flow. This 
corresponds to an expansion about the mean flow and to 
first order yields the following evolutionary equations in 
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where p is the mean density, F r is the mass fraction of 
species r, pui is the momentum in the i th direction, P 
is the mean gas pressure, and e is the internal energy 
per unit mass which, unlike in DT06 and SB08, includes 
both the thermal and turbulent kinetic energy compo- 
nents. As discussed in Scannapieo & Briiggen (2010), 
this formulation allows us to follow the model into the 
highly supersonic regime in which most of the internal 
energy is turbulent rather than thermal. Turbulence af- 
fects the mean flow through the turbulent stress tensor 
Tij and the turbulent viscosity pr which is scaled in the 
energy equation by N £ = 1 and in the mass fraction equa- 
tion by Nf = 1. Finally, the Lagrangian time derivative 
is defined as 

D_ = d_ ~d 
Dt dt 



^ U ^ 8Xn ' 



(6) 



where there is an implied summation over all dimensions. 

These equations depend on the evolution of the eddy 
scale L and the turbulent kinetic energy K. Equations 
that include the diffusion, production, and compression 
of these quantities are 
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(8) 



In eqn. ([JJ the first term represents the diffusion of the 
eddy length scale as scaled by the turbulent viscosity 
Pt and scale factor Nl = 0.5. The second term is the 
primary production term and is proportional to V = 
\/2K and independent of the flow. The third term is the 
growth of eddies due to the expansion and compression 
of the mean flow. Finally, Cc = 1/3 is a constant in the 
model and is determined by mass conservation in eddies 
as they are compressed. In eqn. (j8j), which parallels eq. 
(j5j), the first term is the diffusion of turbulent kinetic 
energy and is scaled by pt/N £ . The second term is the 
work associated with the turbulent stress which drives 
the KH instabilities. Finally, the third term is a source 
term that drives RT and RM instabilities. 

Note that we assume that Nf ~ N e . Pan & Scanna- 
pieco (2010) studied the efficiency of mixing over a large 
range of high Mach number turbulent flows. By com- 
paring the scalar (e.g. mass fraction) dissipation time 
scale to the time scale for the total kinetic energy loss, 
they find that this ratio does not deviate much from one, 
which validates this choice. 
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TABLE 1 



Coefficient 


Value 


Effect 


Mr-, 


i n 
1 .u 


Diffusion of Species 


N e 


1.0 


Diffusion of Internal Energy 


N L 


0.5 


Diffusion of L 


Cc 


1/3 


Compression Effects 


C B 


0.84 


Buoyancy-driven turbulence 


Cd 


1.25 


Drag term on K 


Ca 


2.0 


Atwood Number 


C P 


2/3 


Turbulent Pressure 




1.0 


Turbulent viscosity 




variable 


KH growth scaling 



The primary source term for RT and RM instabilities 
is Sk-, which is represented by and defined by DT06, 



Sk = pV 



paxi 



V 2 



(9) 



where the coefficients Cp = 0.84 and Cp = 1-25 are 
fit to turbulence experiments. Physically, turbulent en- 
trainment is described by Cb which reduces any density 
contrasts, and Cp is a drag coefficient that describes the 
dissipation of turbulent energy when the average length 
scale is proportional to L. Likewise, V = V2K is the 
average turbulent velocity, P is the pressure, p is the 
density , and Ai describes the Atwood number in the 



nth 



direction. This is determined by, 



A- = 



P+ 



C A 



dp 



p + L\dp/dxi\ dxi 



(10) 



where Ca = 2 is a constant of the model, p + and p- are 
the densities on the front and rear boundaries of a cell 
in the ^-direction. 

Additionally, and unlike DT06 and SB08, we include 
the full Reynolds stress tensor, constructed from mean 
velocities: 



nj = CpS[jpK - p T t K h 



dui duj 2^ duk\ 
dxj dxi 3 1J dxk J 



where 5-^ is the Kronecker delta function, pp is the tur- 
bulent viscosity, tkh is a function of the local Mach num- 
ber which is calibrated to produce the correct KH growth 
rate as discussed below, and Cp is a constant. The first 
term is the isotropic turbulent pressure and has a trace of 
2pK when Cp = 2/3. The second term is the deviatoric 
tensor which has a zero trace (note the implied summa- 
tion over all dimensions) and is the primary source of 
shear instabilities. 
Finally, the turbulent viscosity is calculated as. 



p T C^pL\ 



(12) 



where C M = 1 is a constant. Table [T] summarizes all 
model coefficients, their values, and their effects. 

The numerical implementation of this model is divided 
into five steps: 

1. Update the velocities using the turbulent viscosity 
in the fourth-order PPM solver before the turbu- 
lence package is called during the hydrodynamic 
step. 



2. Calculate the Reynolds stress tensor and update 
the turbulent kinetic energy, K, as in eqn. (jSJ). 

3. Use the updated value for K and actualize the dif- 
fusive mixing terms in eqns. (j3]), (j5j), (0, and (j8]). 

4. Compute the contributions from the source terms 
as: 

(a) Calculate V = V2K. 

(b) Add the pV term to the L equation and use 
a leapfrog approach to add the source term 
(Sk/ pV) to the V equation. 

(c) Write K back asK = V 2 /2 a'nd update the 
turbulent viscosity as in eqn. ([T2]) . 

5. Enforce a minimum time-step from turbulence as 
dt < (A 2 //i T )/6, where A is the minimum between 
dx, dy and dz in a given cell and pp is the turbulent 
viscosity in that cell. The minimum per block is 
calculated where a block in FLASH is composed of 
nx x ny x nz cells. Finally the global minimum is 
calculated across all blocks. 

2.2. Sub- Grid Turbulence Model Tests 
2.2.1. Ray leigh- Taylor Shock Tube Test 

To verify the implementation of our model, we recon- 
structed the RT problem as described in §5 of DT06 (and 
§3.1 of SB08). Initially a 1 cm region was filled with two 
7 = 5/3 fluids with constant density, p\ = 1.0 gem -3 in 
the region from x = 0.0 to 0.5 cm and p2 = 0.9 gem -3 
frormc = 0.5 to 1.0 cm. A gravitational field acted in the x 
direction with a constant acceleration of 9.8 x 10 8 cm s -2 
or 10 6 times the Earth's gravity. The initial tempera- 
ture profile was calculated so that both fluids were in 
hydrostatic equilibrium and so that at x = 0.5 cm the 
temperature of the lower density fluid was T2 = 50 K 
and the temperature in the higher density fluid was T\ 
= 45 K. Finally, to test the mass fraction diffusion, we 
initialized each side with different generic mass fractions 
with atomic masses equal to hydrogen. 

Despite being described as a 1-dimensional problem, 
to test our implementation in FLASH we set up a 2- 
dimensional 50 "block" by 1 "block" region. Each block 
represented 8x8 simulation cells. Each test was allowed 
to refine up to a l re f = 4 based on density and pressure 
profiles. This led to an initial cell size of 1.0/50/8/2 3 = 
3.1 x 10 -4 cm at the interface and a minimum resolution 
of 2.5 x 10 -3 cm. 

Although there is an explicit turbulent time step that 
must be enforced, this implementation works very well 
with the AMR hydrodynamic time step imposed in 
FLASH. Initially when the center is fully refined, the 
hydrodynamic time step is shorter than the one imposed 
by turbulence. As the turbulent viscosity grows the time 
steps become comparable, however, because of the dif- 
fusion associated with turbulence the density contrast is 
smoothed. This allows the AMR to derefine these areas, 
which in turn increases the turbulent time step. Finally 
the turbulent time step reaches an equilibrium of 1 /6 of 
the hydrodynamic time step at the end of the simulation 
after the whole volume has reached the lowest refinement 
level. 
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As described in DT06, this problem has an analytic 
solution for the evolution of K and L, 



K(x,t)f*K (t)[l 



x 2 \ 
h?(t)) 



L(x,t)^L (t)Jl 



h 2 (t)> 



(13) 
(14) 



where h(t) is the scale length for the interpenetrating 
fluid and is defined as, 



h(t) = a b A(0)gt 2 , 



(15) 



here = 0.06, A(0) is the initial Atwood number = 0.05, 
g is the gravitational acceleration and is set to 9.8 x 10 8 
cms -2 , t is time in seconds, Ko(t) = (dh/dt) 2 /2, and 
L (t) = h(t) 2 /2. Initially both K and L were set to a 
small values throughout the simulation except near the 
interface where we set both K and L to the analytic 
values at time of 50 /as. 

We find that our model matches the expected profiles 
for K and L at all times and as expected L and K evolve 
as oc t 2 . As K and L quickly increase, the mixing layer 
also increases which promotes rapid mixing between the 
two fluids which can be seen in the species profiles. We 
carried out the same test expect with a sharper density 
contrast with p 2 = 0.8 gm/cm 3 and T 2 = 40 K. The 
results are shows in the second column of Fig. [TJ Again, 
there is excellent agreement with the expected analytic 
profiles. 

2.2.2. Shear Flow Test 

To test the ability of our model to accurately model 
subsonic shear flow mixing we adapted the shear flow test 
problem described in §3 of Chiravalle (2006). We began 
the problem with an initial velocity shear discontinuity 
at the origin. Left of the origin we set the y- velocity to 
7.8 x 10 4 cms" 1 (Mi = 0.46) while on the right we set 
the velocity to 1.09 x 10 5 cms" 1 (M 2 = 0.66). Pressure 
and density were held constant across the full domain at 
1.72 x 10 10 erg cm -3 and 1.0 g cm -3 respectively. To 
study mass fraction diffusion, we again initialized each 
side with different mass scalars with identical properties. 

In this case, the shear layer is expected to grow linearly 
with time as 

5 = 0.181 Avt, (16) 

where S is the width of the mixing layer, Av is the 
difference in velocity across the interface, and t is the 
time (Chiravalle 2006). Unlike the test in Chiravalle 
(2006), we also added a small initial shear layer of size 
^init = 0-1 cm split equally through the interface with 
K = 0.02 (Av) 2 and L = 0.2 tf init . With this setup, we 
were able to vary tkh as a free parameter to approxi- 
mate the expected growth rate. After 150 /is, the ex- 
pected width of the shear layer is 0.923 cm, we find that 
tkh = 0.20 reproduces this result, as shown as the top 
line in Fig. O This value is close to the one suggested by 
DT06 (r KH -0.1). 

2.2.3. Supersonic Shear Test 



To extend our model into the supersonic regime, we 
compared our mixing layer widths to those obtained ex- 
perimentally by Papamoschou & Roshko (1988), who 
measured the growth rate of the shear mixing layer as 
a function of Mach number by forcing two fluids across 
each other at different relative speeds. Defining the con- 
vective Mach number as 



M c 



\Ui-U2\ 
a\ + a 2 



(17) 



where a\ and a 2 are the sound speeds and U\ and Z7 2 
are the fluid velocities in region 1 and 2 respectively, Pa- 
pamoschou & Roshko (1988) found that as M c \ increases, 
the mixing layer quickly decays asymptotically to 20% of 
the subsonic mixing layer width. To match this behavior 
at high Mach numbers we allow our variable tkh to vary 
depending on a 'local' Mach number which we define as 



Mi 



l(VxG)|L 



(18) 



where |V x u| is the magnitude of the curl of the mean 
velocity field, L is the local eddy scale, and c s is the local 
sound speed. This local Mach number approximates M c \ 
and we use it to scale tkh as 

(0.2 Mi < 0.3, 

tkh = < 0.2 - 0.65(Mi - 0.3) 0.3 < M Y < 0.6, (19) 
[o.00575 0.6 < Mi. 

To compare the result of this approach to the exper- 
imental measurements of S as a function of time, we 
adapted the subsonic shear test from §3.2 by changing 
the velocity on one side of the interface to match the 
the expected M c i. As above, we also initialized a small 
shear layer of size Si n a = 0.1 cm split equally through 
the interface with K = 0.02 (Av) 2 and L = 0.2 tf init . K 
and L are initialized to zero everywhere else. Finally, by 
changing the initial Av (and thus M c i) and evolving for 
the appropriate time we measure the width of the final 
mixing zone. 

Fig. [2] shows the results of selecting a variety of values 
for M c i and varying tkh- Each set of points in this figure 
gives the mixing width measured as the distance between 
the two points which correspond to 1% and 99% of the 
velocity difference. Note that this is discretized due to 
this definition and the spatially discrete (AMR) structure 
of FLASH. Also in Fig.[2]rhe solid lines are the theoretical 
mixing widths of the form 



S = Sq + 0.181 Avtk, 



(20) 



were So is the initial mixing width, Av is the velocity 
difference between the two fluids, t is time, and A: is a 
constant between and 1. 

In this figure the simulation time has been normalized 
by the evolution time, defined as 0.812 cm/Av, the time 
required for each case to reach a final mixing width of S 
= 0.923 cm if k were equal to 1. If tkh is scaled correctly, 
then when fitting the mixing width using eqn. ([20]) k will 
equal the expected mixing layer widths divided by the ex- 
pected width as a function of Mach numbers as given by 
Papamoschou & Roshko (1988) in their Fig. 16. As Fig. 
[2] shows, our model reproduces the expect linear growth 
extremely well across a range of Mach numbers. This 
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Fig. 1. — Evolution of the shock tube tests. The p2 = 0.9 g cm -3 case is given in the first column and the p2 = 0.8 g cm -3 case is given 
in the second column. Top left panel: Profiles of turbulent length scale at 50 (red), 100 (blue), 200 (green), and 300 (magenta) fis. In 
each case, the dotted lines are the analytic solution and the solid lines are the simulation results. Second left panel: Profiles of the kinetic 
turbulent energy of the same case at the same times as the top panel. Third left panel: Density profiles at the same times. Fourth left 
panel: Temperature profiles. Bottom left panel: Profiles of species mass fractions. The dotted lines show the mass fraction of the species 
that was initially on the left side and the solid for the species on the right side. Right panels: same as the left except the profiles correspond 
to 50 (red), 100 (blue), 150 (green), and 200 (magenta) ps respectively. The x-axis and y-axis scales are the same in both columns. 



TABLE 2 

Parameters used for the Fig. [2] The first 2 columns are 
the initial velocities on either side of the interface in 
units of mach numbers, the 3rd column is the convective 
Mach number, the 4th column is the evolution time for 
each model, the 5th column is the final mixing width, and 
the 6th column is the parameter k as defined in eqn. [20j 



Mi 


M 2 




r (/zs) 


6 (cm) 


k 


0.46 


0.66 


0.10 


150 


0.923 


1.0 


0.46 


1.46 


0.50 


30 


0.48 


0.52 


0.46 


1.96 


0.75 


20 


0.32 


0.32 


0.46 


3.46 


1.50 


10 


0.15 


0.16 



range is much larger that seen our simulations, which 
have typical Mach numbers between 0.3-0.7. Table [2] 
summarizes the initial setup parameters for the results 
in Fig. [2] as well as the final mixing widths and values for 
k. 

The growth rate of a shear layer is dependent primar- 
ily on the velocity and density ratio on either side of 
the shear discontinuity. This dependence has been stud- 
ied by numerous authors (e.g. Brown & Roshko (1974), 
Slessor et al (2000) and references within) and as shown 
by Soteriou & Ghoniem (1995) is small and for a given 
velocity ratio does not alter the growth rate very much. 



2.3. Radiative Cooling 

Above 10 4 K the cooling function is primarily con- 
trolled by atomic radiation and, at very high tempera- 
tures, bremsstrahlung radiation. These contributions are 
calculated from a table lookup using values calculated 
using CLOUDY (Ferland et al 1998). Here we assume 
Case B recombination and consider only collisional ion- 
ization by use of the "coronal equilibrium" command of 
a metal free gas. Below 10 4 K the cooling is dominated 
by molecular line cooling and metal-line cooling. 

Molecular cooling is described in Paper I, but now in 
addition, we have included metal line radiative cooling 
in the optically thin limit. To simulate metals in our 
simulations, we define a generic mass scalar in FLASH, 
which is advected with any flows. For the cooling rates 
we use the tabulated results from Weirsma et al (2008), 
which assume local thermodynamic equilibrium, and in 
all cases we use standard solar abundance ratios. The ra- 
diative rates are defined over a large temperature range, 
from 10 2 K through 10 9 K. The specific cooling rate for 
a given temperature is found using a table lookup from 
a data file and scaling by the local metallicity. 

3. MODELING OUTFLOW-MINIHALO INTERACTIONS 

Having described the new physical processes, we return 
our attention to the model developed in Paper I. Again 
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Fig. 2. — Expected growth rate at different Mach numbers. The 
y-axis is the width of the mixing layer in cm and the x-axis is the 
normalized evolution time (the simulation time divided by the total 
evolution time 0.812cm/Ai;,). The red lines are the expect mixing 
widths from the k values given in Papamoscho Sz Roshko (1988) 
while the blue points are the measured widths from our model. 
Table [2] summarizes the parameters used for the fits and for the 
model. The range of Mach numbers studied here is much larger 
than the range of Mach numbers found in our simulations, which 
vary between 0.3-0.7. 

we assume a ACDM cosmology with h = 0.7, = 0-3, 
^ A = 0.7, and Q b = 0.045 (e.g., Spergel et al 2007), where 
h is the Hubble constant in units of 100 km s _1 Mpc -1 , 
and f^o, £^a, and are the total matter, vacuum, and 
baryonic densities in units of the critical density. The 
critical density for our choice of h is p cr i t = 9.2 x 10 -30 
g cm -3 . 

We begin with a neutral primordial minihalo, which is 
composed of 24% helium and 76% hydrogen with a total 
mass of both dark and baryonic matter of M c = 3.0 x 
10 6 M . The initial minihalo has a total radial density 
profile given by Navarro et al (1997): 



p(R) 



QqPc 



cx(l + cx) 2 SF(c) 



g cm 



(21) 



where c is the halo concentration factor, x = R/R c , 
R c {= 0.393) kpc is the virial radius, F(t) = ln(l + t) 
- yq^-, and p c = 6.54 x 10 -25 g cm -3 is the mean cluster 
density. The baryonic matter is taken to be in hydro- 
static equilibrium and follows an isothermal radial profile 
with a virial temperature of T = 1650 K: 



Pgas(R) = Poe 



g cm 



(22) 



4 sc (xR viT ) = 2v 2 c [F{cx 
2.16 x 10 



9^ — *3 



where the escape velocity is v 2 
cx(l + cx)~ 1 ][xF(c)]~ 1 and po 

Gravity is treated using the multigrid Poisson solver for 
self gravity of the gas (Ricker 2008) as well an additional 
component of gravitational acceleration due to dark mat- 



ter. The initial metallicity of the halo and surrounding 
gas is set to zero, and the initial values of K and L are 
set to 1% of the total internal energy and one parsec re- 
spectively. All other parameters are left at their fiducial 
values including the background UV radiation field (J21 
= 0.0). 

The outflow is approximated by the Sedov- Taylor blast 
wave solution. We assume that the minihalo is at a dis- 
tance R s = 3.6 kpc and that the shock has a velocity of 
v s = 225 km s _1 . By the time the shock reaches the mini- 
halo it will have entrained a total mass of M S)to tai = 4.4 
x 10 T M with an associated surface density of a s = 2.6 
x 10 5 M kpc -2 . The input energy for this outflow was 
derived from SNe from the host galaxy, and was taken to 
be E = (eEss) erg where e is the wind efficiency which 
is obtained from the fraction of the SNe energy funneled 
into the outflow and is set at our fiducial value of e = 0.3 
and E55 is the total SNe energy in units of 10 55 erg. 

We initialized the left boundary with the same initial 
properties as our previous models. As in the rest of the 
simulation domain, we set K and L to 1% of the inter- 
nal energy and one parsec respectively. To determine the 
initial metallicity of the shock, we followed the analysis 
from Scannapieco et al (2004). We infer that roughly 
2 M worth of metals are produced per 10 51 ergs of en- 
ergy in both Type II and pair instability supernova from 
Population III stars (Woosley & Weaver 1995; Heger & 
Woosley 2002). If we assume that half of these met- 
als escape from the host galaxy and are funneled into 
the outflow, we can expect a total mass of metals of 
^metai = 10 4 E55 M . This leads to a metal fraction of 
the shock of X meta i = M meta i A4 5 totai = 0.12 Z , which 
we use as our fiducial value. 

4. RESULTS 

Our simulations were carried out in a rectangular box 
with an effective volume of 3.2 x 10 9 pc 3 . The y- axis and 
z-axis were both 1170 pc and ranged from (-585,585) pc. 
The x-axis was twice as long, 2340 pc, and ranged from 
(-585,1170) pc. The minihalo was centered at (0,0,0) pc 
and the shock originated from the left boundary with a 
velocity along the positive x direction. Both density and 
pressure were used as the refinement /derefinement vari- 
ables, and we also forced derefinement after 7 Myrs in re- 
gions with density less than 3.26 x 10 -26 g cm -3 outside 
a central sphere of radius 324 pc centered at (0,0,0) kpc. 
This had the advantage of derefining unimportant blocks, 
which was especially important in the runs without sub- 
grid turbulence, as turbulent mixing tends to smooth the 
density gradients, allowing the AMR to naturally dere- 
fine. 

A summary of the runs performed are given in Table 
[3j We label them by whether they were high, medium, 
or low resolution (H, M, or L) and whether they used the 
subgrid turbulence package (WT or NT). In our new sim- 
ulations metal cooling was always included, and to asses 
the impact of this cooling we also include the fiducial run 
from Paper I (HBN), noted by the asterisk, which is used 
to compare with run HNT. 

4.1. Hydrodynamic Evolution 

Figures [3] and H] show the evolution of the minihalo from 
initial setup through the shock interaction and to the 
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Fig. 3. — Evolution of run HWT from t = to the time at which the shock completely surrounds the cloud. Each image shows an 
x-y slice through the center (z=0) of our simulation volume. The first column shows logarithmic density contours from 10 - g cm -3 to 
10 -21 g cm -3 , which correspond to number densities between n « 10 -2 cm -3 and 10 2 cm -3 . The second column shows the logarithmic 
temperature contours from 10 K to 10 8 K, the third column shows the logarithmic H2 mass fraction contours between Xh 2 = 10 -8 to 
10 _1 , and finally the fourth column shows the logarithmic metal mass fraction contours between Z = 10 -4 Zq to 10 -0-5 Zq. 
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Fig. 4. — Evolution of run HWT from the time at which the shocks meet at the back of the cloud (t = 6.797 Myrs), to the time at which 
the reverse shock passes through the cloud (t = 7.67 Myrs), to the collapse of the cloud (t = 11.9 Myrs), and the end of the simulation 
(t = 14.65 Myrs). Columns and rows are the same as in Fig. [3] For this figure we have cropped the individual images along the x and y 
axes for clarity. Until t= 7.67 Myrs the molecule and metal distributions closely follow each other, but at later times molecule formation 
is enhanced near the core of the cloud due to the reverse shock, which does not carry metals. 
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TABLE 3 

Summary of the simulations. 



Name Z re f Resolution (pc) Turbulence Metal Cooling 



LWT 4 

LNT 4 

MWT 5 

MNT 5 

HWT 6 

HNT 6 

HBN* 6 



18.22 

18.22 

9.11 

9.11 

4.55 

4.55 

4.55 



Y 
N 
Y 
N 
Y 
N 
N 



Y 
Y 
Y 
Y 
Y 
Y 
N 



final collapse of the cloud, focusing on several important 
stages of evolution throughout these figures. The first 
row in Fig. [3] shows the initial minihalo. Because the gas 
consists of neutral hydrogen and helium, it is unable to 
cool on its own. Instead it remains in hydrostatic balance 
with a free fall time of 



is 



3tt 
32Gp 



100 Myrs, 



(23) 



that is roughly equal to the sound crossing time. 

As the shock interacts with the minihalo, it begins to 
heat and ionize the neutral gas. The ionized gas then re- 
combines and starts to catalyze the formation of H2 and 
HD. This results in a 'hollow' distribution of molecules 
that is concentrated in the interacting regions at the 
front and sides of the cloud. The metal distribution 
closely follows this molecular distribution, because the 
shock is able to move faster in the less dense portions of 
the cloud. The shock fully envelopes the minihalo in a 
'cloud-crossing' time defined by Klein et al (1994) as 



tc. 



2Rc 



• 4.6 Myrs, 



(24) 



which occurs at about t « 6.5 Myrs after the beginning 
of our simulations and is shown in the last row of Fig. [3] 
The time scale for H2 formation is given by Glover et 
al (2008) as 

«* = itn <25) 

where Xh 2 and X e are the mass fractions of H2 and elec- 
trons respectively, n is the number density, and k\ is the 
reaction rate for the formation of H - a key react ant in 
formation of H2 . When the shock first interacts with the 
cloud, this time scale is very short since n ~ 1 cm -3 and 
the fraction of electrons (X e ) is relatively large. However 
the cloud quickly reaches an abundance of Xh 2 ~ 10 -4 
as shown in the third column of Fig. 03 This fraction con- 
tinues to grow as the cloud is surrounded, but does so at 
a much slower rate as Xh 2 increases and X e decreases. 
Again, the distributions of molecules and metals closely 
track each other. 

At « 7 Myrs, after the shock meets on the back of the 
cloud, it creates a reverse shock, which begins to catalyze 
molecule formation, as shown in the top row in Fig. HJ It 
is here that the metal and molecule distribution diverge 
as the metals have not had the time to diffuse in the 
interior of the cloud. We define a turbulent mixing time 
scale as 

- £ 

tmix " (p T /p) S ' ( 6) 



where d is the distance over which the metals are diffused, 
\±t is the turbulent viscosity, and p the local density. By 
looking at the third and fourth columns of Fig. |4] and 
comparing the distributions of H2 and metals at t = 7.67 
Myrs, it is obvious that metals are deficient along the x- 
axis at y ~ 0.3 kpc. If we approximate the distance that 
metals need to diffuse through as d ~ 10 pc and estimate 
the turbulent viscosity around the collapsing cloud to its 
post-shock value pr / P ~ 25 pc 2 Myr -1 then the mixing 
time scale is ~ 4 Myrs. Thus at « 7 + 4 = 11 Myrs, 
the distributions of molecules (shown by H2) and metals 
once again follow each other, as can be seen in the third 
row of Fig. [U 

The fourth row of Fig. [H shows the final state of the 
cloud. Most of the material is found in a small dense rib- 
bon that is stretched along the x-axis and extends many 
times the initial virial radius away from the center of 
the dark matter halo. This material is now much colder 
than it started with typical temperature around 100 K. 
The H2 mass fraction abundance of this ribbon is around 
10 -2,2 , and it has a metallicity of about 10 ~ 2 Z . 

4.2. Model Dependencies 
4.2.1. Effect of Turbulence 

Turbulence has two primary effects in our simulations: 
the diffusion of metals from the shock into the minihalo 
and the smoothing of sharp density contrasts. Fig. [5] 
compares the difference between runs with (HWT) and 
without (HNT) our subgrid model for turbulence. As 
expected many of the sharp density features found in 
the model without turbulence are diffused away in the 
model with subgrid turbulence. This is seen in the late- 
time panels in Rows 1 and 2 in Fig. [5] where there was 
a prominent lower density feature in HNT (at x ~ 0.65 
kpc) that is not seen at all in HWT. Also absent are the 
smaller density features at the far end of the simulation 
domain along the x-axis. Although, in both cases the 
general result is the same: much of the mass has formed 
into a dense ribbon along the x-axis. 

The most striking difference between the two simula- 
tions is the metal distribution. HWT shows a very uni- 
form metal abundance in the final cloud. Almost every 
portion has a final abundance of Z w 10 -2 Z except 
for a low density region near the initial center of the 
cloud. HNT however shows a much more uneven distri- 
bution with two regions of low metallicity; one at 0.15 
kpc and the other at 0.6 kpc. However, it is important 
to note that in both cases (HWT and HNT) the densest 
regions in both models have essentially the same final 
metal abundance, and thus the spread in metal distribu- 
tion in the stars that are formed would be small with or 
without the inclusion of subgrid turbulence. 

4.2.2. Effect of Metal- Line Cooling 

At temperatures below T w 10 4 K the primary 
coolants are molecules and low-energy metal-lines. 
Therefore the total cooling is expected to be strongly 
dependent on differences between metal and molecule 
abundances. Figure [6] shows the difference between a 
fiducial model with (HNT) and without (HBN) metal 
cooling, where HBN is taken from Paper I. There is little 
difference between these runs aside from slight changes 
in the positions of small structures. In both cases the 
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Fig. 5. — Late time comparison between a run with subgrid turbulence (HWT) and a run without it (HNT). Each column represents a 
different snapshot in time. The top two rows shows logarithmic density contours from 10 -26 to 10 -21 g cm -3 . Rows 3 and 4 show the 
logarithmic contours of H2 mass fraction from 10 -8 to 10 -1 . Finally, Rows 5-6 show the contours of metallicity in units of solar metallicity 
between 1CT 4 - to 1CT - 5 Z©. In each set of rows the model with sub-grid turbulence (HWT) is on top of the model without it (HNT). 
Each image is a slice through the center of the domain along the z-axis. 



same dense knots are found in essentially the same places. 
Furthermore, the abundances of molecular coolants are 
essentially identical and are not affected by the inclu- 
sion of metals. This suggests that metal cooling is not 
as important as molecular cooling, because otherwise the 
abundance of molecules in run HNT would be lower than 
that in run HBN due to the temperature dependence in 
the molecule formation rates. 

The time scale for H2 cooling can be estimated as 



m 2 



l.bnkT 



n e n H2 A(T) H ,H 2 



(27) 



where n is the total number density, k is Boltzmann's 
constant, T is the temperature, n e and nn 2 are the 
number densities of electrons and H2 respectively, and 
A(T)h,h 2 is the cooling rate as a function of tempera- 



ture. Similarly for metal cooling the time scale is 

l.bnkT 



s, 



(28) 



where n e and nn are the number densities for electron 
and hydrogen and A(T)m is the cooling rate for metals. 
The ratio of these time scales is then 



TH* = X e A H2 A(T) M Z 

tm Xft 2 A e A(T)h,h 2 Z 



(29) 



where the number densities have been replaced using X{ 
= TiiAij pN a- To get an idea of the relative importance 
of each cooling method, we take representative values for 
these variables at t = 7.67 Myrs, an important point in 
the evolution of the cloud as it begins to collapse. We find 
that th 2 /tm ~ 10 -2 which means that at this moment 
molecular line cooling is more important than metal line 
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Fig. 6. — Comparison between a run with metal cooling (HNT) and a run without it (HBN). Rows 1 and 2 show the logarithmic contours 
of density for HNT and HBN while Rows 3 and 4 show the logarithmic contours of H2 mass fraction. Apart from the slight differences in 
the positions of similar features there is very little difference between the runs. 



cooling. Overall, the inclusion of meal line cooling does 
not make a significant difference in the final state of the 
cloud. 

4.2.3. Effect of Resolution 

Finally, we study the impact of resolution on our re- 
sults. Fig. [7] shows the result of our models over a range 
of maximum resolution levels. The left two columns show 
the results for models with turbulence while the right two 
columns show the results without turbulence. Generally, 
the outcome is independent of resolution. In the case of 
density the higher the resolution the smaller and more 
dense the final cloud becomes. However, the structure 
of the cloud is the same: it has been stretched into a 
ribbon and moved out of the dark matter halo. Further- 
more, the metallicity distribution is almost indistinguish- 
able between the runs with different resolution levels. For 
each choice of resolution level the final cloud is enriched 
to w 1(T 2 Z/Z , both for runs with and without subgrid 
turbulence. 

Not shown is the difference in molecular abundances. 
Here there is a resolution dependence on the formation 
time scale for molecular coolants, as discussed in Paper 
I. However, the amount formed is always sufficient to 
cool the cloud enough to produce the same final out- 
come. The final abundances are nearly identical at each 
refinement level except at the lowest resolution, which 
is different only outside of the dense regions in the final 
ribbon. Aside from features that do not affect the final 
distribution of star-forming gas, our results are indepen- 
dent of resolution from / max = 4 to / max = 6. 

4.3. Stellar Clusters 



The final distribution of the dense clumps evolves over 
100s of Myrs, a much longer time scale than the shock- 
minihalo interaction itself. To evolve our simulation 
over such a long time scale, we adopt the simple one- 
dimensional procedure described in Paper I. Here we sub- 
divide the x-axis into 100 evenly spaced bins from x = 
kpc to x = 1.4 kpc. The mass for each bin is calculated 
by summing up the gas density from the simulation in a 
cylindrical volume with a radius of 24 pc and length of 
the bin. We also calculate the velocity of each new bin 
by summing together the momentum from each cell that 
goes into it and dividing by the total mass. 

This distribution is then evolved using a simple nu- 
merical model in which we assume all motion is along 
the x-axis and that pressure is not important at late 
times. The acceleration of each point is calculated by 
adding together the gravitational acceleration from all 
other particles as well as the gravitational acceleration 
from the dark matter halo. A leapfrog method is used to 
calculate the updated position and velocity from the cal- 
culated acceleration and updated velocity respectively. 
Finally, if a bin is evolved past the one in front of it, 
we merge these two by summing their mass and using 
conservation of momentum to compute its velocity. 

The model is evolved for 200 Myrs past the end of the 
simulation, and the results are shown in Fig. [8] for runs 
HWT and MWT. As the gas continues to move outward 
the particles begin to attract each other and merge. By 
50 Myrs after the FLASH simulation, most of the parti- 
cles have merged and the motion of the remaining parti- 
cles is purely ballistic. This can be seen in the top and 
middle panels of Fig. [8] as the lines begin to overlap each 
other. 

At the final time of 200 Myrs after the end of the sim- 
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Fig. 7. — The impact of maximum resolution on our runs with (left two columns) and without (right two columns) subgrid turbulence. 
The top three rows show the logarithmic density contours at t = 7.67 Myrs and t = 14.65 Myrs and the bottom three rows show the 
logarithmic metallicity contours in units of solar metallicity. Rows 1 and 4 show the respective contours at the highest resolution with l max 
= 6 (runs HWT and HNT), Rows 2 and 5 at £ ma x = 5 (runs MWT and MNT), and Rows 3 and 6 at Z ma x = 4 (runs LWT and LNT). 



ulation, we can identify one primary clump with a mass 
over 6.0 xlO 4 M with two smaller clumps with masses 
of 4.0xl0 4 M and 2.5xl0 4 M respectively. The two 
largest masses are far outside the dark matter halo, and 
are no longer bound to it. 

Comparing the solid (HWT) and dotted red line 
(MWT) we conclude that regardless of the resolution 
level the same general conclusions hold. There are some 
differences between the two resolution levels, primarily 
in the speed and position of the final clumps formed, but 
these are minor. In each case a nearly identical distribu- 
tion of stellar clusters is formed: one large cluster with 
a few smaller neighbors. 

4.4. Relation to Halo Globular Clusters 

As discussed in Paper I, shock-minihalo interactions 
could be a source of halo globular clusters, as the knots 
will continue to collapse into dense stellar clusters and 
the longest-lived of these stars will survive to the present 
day. In fact the clusters in our simulations are very dense 
(n ~ 10 2 cm -3 ) and expected to become gravitationally 
bound to larger structures that form over cosmological 
time. 

There are several other properties of our simulated 



clusters that support this connection. Globular cluster 
masses are well defined by a Gaussian in logio(MGc /Mq) 
with a mean mass of 10 5 M and a dispersion of 0.5. This 
represents a spread in globular cluster masses of 3.0 x 10 4 
M to 3.0 x 10 5 Mq, which spans the range of the stel- 
lar cluster formed in our simulations. However the final 
cluster masses may depend on the initial mass in the 
minihalo a parameter study is required to determine this 
dependence. This is underway and will be presented in 
a future paper. 

The lower mass limit for globular clusters seems to 
be set by several destruction mechanisms, which include 
mechanical evaporation (eg., Spitzer & Thuan 1972) and 
shocking as the cluster passes through the host galaxy 
(eg., Ostriker et al 1972). On the other hand, the upper 
mass limit seems to be a property of the initial popu- 
lation (eg., Fall & Rees 1985; Peng & Weisheit 1991; 
Elmegreen 2010). In fact, the upper limit of « 10 6 M 
roughly coincides with the virial temperature of T^IO 4 
K which corresponds to the limit at which atomic cooling 
becomes inefficient. 

The metallicity of halo globular clusters provides an- 
other constraint. The intracluster metallicity distribu- 
tion is well described by a Gaussian with a mean value 
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Fig. 8. — Long-term evolution of the distribution of particles after 
the end of the FLASH simulations. The x-axis gives the cumulative 
mass in units of a solar mass. The top panel shows the mass of 
each particle, the middle panel shows the velocity of each particle 
in units of kilometers per second, and the bottom panel shows 
the position of each particle in units of kpc. The solid lines are 
taken from run HWT, and the green line shows the initial profile 
at tf = 14.65 Myrs, the blue line shows the profile at 50 Myrs 
later, cyan shows the profile at tf + 100 Myrs, magenta shows 
the profile at tf + 150 Myrs, and the red line shows the ty+200 
Myrs. By tf + 50 Myrs most of the merging has finished and the 
largest clumps have formed. The dotted red line shows the tf + 
200 Myrs distribution for run MWT and illustrates the difference 
between the same model at different resolutions. Although some 
minor differences are present, the same conclusions can be drawn 
from both runs. 

of [jf] ~ -1.6 and a dispersion of 0.3 (Zinn 1985; Ash- 
man & Bird 1993). The metallicity dispersion within 
a given globular cluster is small, usually within 0.1 dex 
(eg., Suntzeff 1993), although in some cases additional 
late-time star-formation from reprocessed material may 
complicate the final observed distribution (e.g. Piotto et 
al 2007; D'Ercole et al 2008; Bekki 2010). Our model re- 
produces the expected mean abundance (Z ~ 10 -2 Z ) 
with an extremely homogeneous distribution given our 
initial abundance in the shock. 

Note that our model keeps track of the velocity, \/2K, 
and eddy turnover scale, L, of bouyancy-driven and 
shear-driven turbulence, and assumes that below these 
lengths scales the flow will behave as fully developed 
turbulence. In this case, as studed in detail in Pan et 
al (2010), the mixing of metals is driven by a cascade 
process similar to that of the velocity field. Using direct 
numerical simulations, Pan et al (2010) showed that over 
a large range of Mach numbers which span the values in 
our in shock minihalo interactions, metals are mixed in 
on a time scale which is close to the time scale for energy 
decay, and that the dependence of this mixing time on 
the length scale at which pollutants are injected is also 
consistent with this cascade picture. 

The final constraint comes from the study of tidally 



disrupted stars in globular clusters. Observations of glob- 
ular clusters have shown tidal tails of stars actively be- 
ing removed from these systems (Irwin & Hatzidimitriou 
1993; Grillmair et al 1995) and have been used to ar- 
gue that globular clusters do not contain dark matter 
(Moore 1996; Conroy et al 2010). If globular clusters 
reside within dark matter halos these tidal stars would 
be suppressed, which suggests that they formed via a 
different mechanism than galaxies. 

Our model naturally reproduces this feature as the 
shock moves almost all of the minihalo gas outside of the 
dark matter halo. As initially inactive gas interacts with 
the shock it is not only enriched by the metals to levels 
expected in globular clusters, the neutral gas is ionized 
and molecular species form via non-equilibrium chemical 
reactions. Not only does the gas have two methods of 
cooling below the T « 10 4 K threshold, but the shock 
imparts enough momentum to remove it from the dark 
matter halo. The result is a population of stellar clusters 
with properties expected of halo globular clusters. 

Furthermore, the clusters formed in our simulations 
should be observable with the James Webb Space Tele- 
scope (JWST) as ribbons of dense star- forming gas that 
are elongated versions of the more compact super-star 
clusters seen at low redshift (e.g. Turner et al 2000; 
Turner & Beck 2004). At 2=10 JWST will have an angu- 
lar resolution of ~ 0.1" and our stretched ribbon of mate- 
rial will span « 0.3", covering nearly three resolution el- 
ements. We can also estimate the apparent magnitude of 
an idealized cluster with a mass of 10 6 M . By assuming 
a single burst of star formation we calculate the appar- 
ent AB magnitudes in the rest frame U, B, and V bands 5 
Myrs after the burst of 32.0, 33.0, and 33.1 respectively 
(Bruzual & Chariot 2003). Given the integration times 
appropriate for deep fields, these clusters can be detected 
by NIRCAM in the 4-6 jim wavelength range (Gardner et 
al 2006). Although smaller clusters were formed for our 
fiducial model here, an upcoming parameter study will 
probe whether similar interactions with larger minihalos 
will form larger clusters with unique elongated morpholo- 
gies. 

5. SUMMARY 

Cosmological minihalos provide the building blocks 
from which larger structures, such as galaxies, can form. 
These structures require interactions from some outside 
source before they can evolve further, as they are filled 
with pristine, metal-free neutral atomic gas that sur- 
rounds the earliest galaxies. Previous work has centered 
on using ionization fronts to induce star formation (Cen 
2001), although follow up work using 3D hydrodynamic 
simulations showed that this would disrupt the cloud 
rather than induce molecule formation (Iliev et al 2005; 
Shapiro et al 2006). On the other hand, galactic outflows 
provide a method of not only inducing molecule forma- 
tion but also providing metals that can be mixed into 
the cloud, perhaps resulting in dense stellar clusters that 
evolve into today's halo globular cluster population. 

However, this scenario is only feasible if these metals 
were able to be efficiently mixed into the minihalo gas. 
While the metallicity difference between globular clusters 
can be quite large, in a given cluster the dispersion of 
(Fe/H) is quite small. It is usually less than 0.1 dex, or 
a factor of ~ 1.25 (Suntzeff et al 1993) in the majority 
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of clusters that do not show late-time star formation, 
which is likely from reprocessed material (e.g. D'Ercole 
et al 2008). We have shown here that in our fiducial 
interaction, most of the primordial gas is enriched to a 
near constant value of Z « 10 ~ 2 Z over a time that 
roughly coincides with the H2 formation time. 

The chemical homogeneity is typically explained by a 
previous generation of supernovae, which either enriched 
the cloud from within, the 'self-enrichment' case (e.g., 
Elmegreen & Efremov 1997; Bromm & Clarke 2002) 
or enriched the gas from which the cloud formed, the 
'pre-enrichment' case (e.g., Brown et al 1995; Cen 2001; 
Nakasato et al 2000; Beasley et al 2003; Li & Burstein 
2003). Both of these scenarios have problems. In the 
'pre-enrichment' case it is unknown what this previous 
generation was, why it did not play a more significant 
role in the evolution of the cloud, and how it can en- 
rich the cloud on such small time scales. However, the 
'self-enrichment' case is even more dire as it requires su- 
pernova exploding inside the cloud without completely 
unbinding the system. In fact, recent 2D simulations of 
such explosions have found that they do indeed destroy 
the minihalo over a wide range of masses and supernova 
energies (Whalen et al 2008b). 

Here we have provided a model that explains many 
of the properties of present-day halo globular clusters, 
through the interactions of a minihalo with a galaxy out- 
flow. In our picture, the outflow performs three primary 
jobs. First, it provides the momentum to move the pri- 
mordial minihalo gas out from its dark matter halo. Sec- 
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